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Abstract 

Numerical solution of nonlinear stiff initial-value problems by a per- 
turbed functional iterative scheme are discussed. The algorithm does not 
fully linearize the system and requires only the diagonal terms of the 
Jacobian. Some examples are presented. 

1. Introduction 

Stiff ode’s (ordinary differential equations) are of special significance 
in chemical kinetics, where decay of one component of the solution could 
happen much faster than other components. For numerical solution of these 
problems, it is generally found that some transient components which are 
negligible in comparison with other components, restrict the step sizes of 
the explicit methods to be very small in order that stability of the numeri- 
cal solution may be maintained. Also, in order to understand the proper 
nature of the solution, the transient effects of the fast decaying terms 
should be revealed by the numerical process. This poses a very difficult 
problem. During the last decade, extensive research has been performed by 
mathematicians, engineers, and chemists to reveal the mechanism of stiff 
systems both mathematically and computationally. 



Most conventional explicit methods such as Euler’s method, Runge-Kutta 
schemes, Adams-Bashforth schemes, etc., require very small step sizes so that 
the algorithm may remain stable. Although some attempts have been made to 
extend the stability properties of explicit methods for special types of 
ode's [ll] by far the most common technique for solving stiff systems numer- 
ically is the use of implicit methods which requires the solution of simul- 
taneous equations. In nonlinear cases Newton-type methods also requires the 
evaluation of a Jacobian, a process that has a high arithmetic operation 
count. These requirements can be costly in both computer time and storage 
requirements. 

Let us now pose the problem, and set up an implicit algorithm for the 
solution. 

^ = f(x,t) (1) 

where x = (x^ X 2 ... f = (fl f 2 ... “ ^o' 

Approximating dx/dt by a backward Euler scheme, we find 

n+1 n n+1 . 

X - X = At f(x , 

or 

x'^'*’^ - x" - At • fCx’^^^, = 0 (2) 

This equation must be solved for x ^ assuming that x’^ the value of x 
at the previous time step is known. Then, we may represent this equation 
as 

F(x"‘*‘^) = 0 (3) 

If X = x’^'*'^ t then it becomes 

F(X) = 0 (4) 

which is a nonlinear system. After we have solved this system for X, we 
update the system by replacing x'^ by the computed values of X and 
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again solve another nonlinear system of the form (4), etc. The process will 
continue until we reach the steady state which is here assumed to exist. 

At present, there are several methods to solve nonlinear systems of the 
form (4). They may be arranged into two classes; 

(i) Methods which require relatively small computer memory storage and 
operational counts per iteration but have slow rates of convergence (e.g.. 
Point Jacobi, Point Gauss-Seidel) . 

(ii) Methods which require relatively large computer memory storage and 
large operational counts per iteration but demonstrate fast rates of conver- 
gence (e.g., Newton’s method). 

In [1] a functional Iterative scheme has been developed to solve non- 
linear systems. It is obtained by perturbing nonlinear Gauss-Seidel itera- 
tions. The perturbation parameters stabilize the algorithm, control the 
mode of convergence, and, in some cases, can greatly speed up the rate of 
convergence. These parameters are essentially corrective factors for the 
iterates. Thus as the iterations .converge, the perturbation parameters 
are all damped out. The present method has been applied to solve several 
stiff nonlinear ODE's. The results seemed to be encouraging in comparison 
with those obtained by •conventional explicit schemes. The primary objec- 
tive in this article is to show the usefulness of this method with regard 
to the numerical solution of stiff ODEs. A brief discussion regarding 
the details of the algorithm and its convergence properties is also 
included. 

2. The Algorithm 

Given a differential equation 

h(x,x,t) = 0 , xCtp) = Xq 

where x = dx/dt and t is the independent variable, if a difference 
approximation is made, it is reduced to a nonlinear system of the form 
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= 0 , 


J Ij 2j 


(5) 


Fj(xi, X2, 


This may be expressed as F(x) = 0, x = (xj X 2 . . . We assume that a 

solution exists. Let us express (5) as: 


X = Gq(x) (6) 

GqI DCR’^ -> D (r’^ = real n-dimensional space). Nonlinear Gauss-Seidel 
iterations to solve (6) may be expressed as: 

x^) (7) 

where x^ = (xi^ X 2 ^ ... x e D; x.^ = value of x, at the kth level 

n J J 

of iteration; G: DxDCR’^xr’'^^D. 

Let us perturb (7) and i^ite: 


= 0)^ + G(x^'*”^, x^) (8) 

where o)^ = o) 2 ^ ... ^ ^ perturbation parameter. 

Ic f 

If the method converges, lim x = x* and G(x*, x*) = x*, which implies: 

k->“ 

Theorem: 1 A necessary condition so that (8) may converge is that for 
some norm 


lim II II =0 (9) 

k->-<» 


In the element form (8) may be expressed as: 


X k+1 = oj.k + 

j J J 


where 


G^k+l,k 


= Gj(: 


k+1 k+i 


k+i 

j-1 


k 



( 10 ) 


( 11 ) 
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Ic 

The algorithm has been derived in details in [1] and Uj is computed by 
using the following equation; 


(0 . = 

J k+r,k 

l-3jGj 


where 


is computed using (11) , 


( 12 ) 



• • • 


jjk+l Qk+l,k 

j-1’ j ’ 


k 

j+1 



and 


1 .L 1 9G 

9 = 

J J 9x 


9G.1 

J 

8xj 


k+i k+i k+i k+i,k k 

> •••» * 




n 


k*4"l k 

Thus in (12), for , (d^ is computed in terms of quantities known 

a priori. The criterion for convergence is, if at some k. 


max 1 03^ I < G 
J 
J 

where e is positive and arbitrarily small. 

3 . Convergence Analysis 

Given a sequence of scalars k = 1> 2, 

D-element iff 


(13) 


a, is called a 
k 


lim a, a^ ... a, =0 (14) 

12k 

It is obvious that if Vk > K, |aj^| 1 a < 1 equation (14) is satisfied. 

Given a sequence of square matrices of the same type with vari- 

able elements. A, is called a D-matrix iff 
’ k 

^ -9 
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Obviously if = A Vk, A^ 
A sufficient condition that 


is a D-matrix iff A 
Aj^ is a D-matrix is 




is a convergent matrix, 
for some norm 


ll\ll i « < 1 


( 16 ) 


The necessary and sufficient condition so that A^ is a D-matrix is 


max 

i»j 


k } k— 1 y ... 1 

a. . 

i.J 


< e 


( 17 ) 


Ic lC“” 1 • • • 1 

where e is positive and arbitrarily small and a.’, » ••• = an element 

J 

of the product matrix 

\ \-i • • • 

An Example ; Let A^ = diag (a,, aj^, bj) 

A2 = diag (a2, 02, b2) , A3 = diag (a3, b3, 03) 


A4 = diag (04, 34, b4) , etc. 

Let lct|-a<lV >N and a 's and b.'s are bounded and chosen arbitrarily 

‘ j j 3 3 

and a finite number of them are greater than 1 in absolute value. Then 


lim A A ... A - 0 ( 18 ) 

n->.co n n-1 1 

Thus each A is a D-matrix, however none of these matrices is a convergent 
3 

matrix . 

Let us consider an iterative scheme 

k+i k+i k. 1, _ 1 9 

X =G(x ,x), k- 1 ,^,... 

Let for a given x* S D, 

G(x^'*‘^ xS - G(x*, X*) = - X*) + Bj^(x*^ - X*) ( 19 ) 

Both A and B, are matrices with variable elements which change as k 
k k 

changes . 
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-I 


G is called a D-mapplng on DxD iff (I - exists 


Vk and (I - Aj^) is a D-matrix. 

The following theorem may now be proved. 

Theorem; 2 In (8) if G is a D-mapping and x = x* G D, then the scheme of 
iteration converges to x* if lim |oo, | = Q. Furthermore, if 

k 00 ^ 

p || (I - < 1, Vk > K, X* is the unique root in D (p(A) = spectral 

radius of a matrix A) . 

Proof ; 

x^"*”^ - X* = 0)^ + G(x^^^, x^) - G(x*, X*) 

= ( 1 ) + - X*) + Bj^(x - X*) ( 20 ) 


Let ” (I “ ^k^~^ ^k ~ from (20) 




x^^^ - X* = ^k^^^ - X*) 


~ ^ \ \-r*‘ \ \+i ••• 


= (E^E 
k 


k -k-1 ••• \+i^ X) %-i 


3 = 1 


O 


+ 2^ \-i ••• "k \-i ••• s <='” - **> 

3=ko+l 

Now lim 1 03^ I =0 implies for some k > k + 1, |ai^| < c (e is positive 
k-~ - 

and arbitrarily small) . Thus each element is arbitrarily small in 

absolute value. Also, Ej^ being a D-matrix 


lim E, E, . . . E =0 
k^« k k-1 1 - 

and lim E, E, . . . E, ., = 0. Hence, - x*| < e which establishes 

k k-1 ko+1 

convergence. 

To prove uniqueness we assume y* G D is another root. Then 
X* - y* = G(x*, y*) - G(y*, y*) 
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This gives 


X* - y* = E^(x* - y*) 

* 

where lim E, = E.. Thus, (I - 1e.|)|x* - y*l - 0. Since p(|E, |)<lYk > K, 

p(|e I) < 1. Hence, by Neumann’s lemma, (I - Ie^I)"^ exists and is non- 
negative. Hence |x* - y*| ^0 which implies x* = y*. 

More discussions on D-matrices and D-mappings are given in [15]. 

4. Some Comparison with Other Methods 

Implicit methods for solving differential equations usually demonstrate 
better stability properties than most explicit algorithms. When a system of 
nonlinear differential equations is reduced to a nonlinear difference system, 
various iterative schemes are available for numerical solution. The simplest 
methods involve functional iterative schemes (Jacobi or Gauss-Seidel itera- 
tions) which have generally very slow rates of convergence. Newton’s method 
has a quadratic rate of convergence which may be expressed as: 

II - X* II = all x^ - x*l|^ 

This immediately implies that even if 0 < a < 1, if || x° - x* || > 1, the 
method could fail . Thus one basic requirement is that the initial guess 
x° should be sufficiently close to the actual root x*. This implies 
II x° - X* II <1. For many initial-value problem this condition could be 
satisfied. However the main disadvantage is that at each k level of 
iteration and for each iterate x^ a unique Jacobian must be computed which 
is not quite practical and for large systems such computations are very 

expensive. 

The present method may be expressed as a combination of nonlinear Gauss- 
Seidel iterations and Lieberstein’s method as follows; 
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k-f^ 







k+l k-t^ 
*1 



k« k+i 

Xj , X2 , 


j = 1, 2, ... n; k = 0, 1, 2, ... 

Thus if rj is the rate of convergence of Gauss-Seidel iteration and r 2 
is the rate of convergence of Lieberstein ' s method, the rate of convergence 
of the present scheme is riT 2 . Also, if at each level Lipschitz condition 
is satisfied, combining both steps it may be seen: 

-X*|i = b2 ||x'‘-x*|| 

where 0 - 3 < 1. Since 3^ goes to zero at a quadratic speed, the method 
converges even if || " x*l| > 1. For various nonlinear systems this global 

convergence property has been verified computationally [14]. 

Furthermore, whereas Newton* s method requires (n^ + n) functionals to 
be computed at each iteration level, only 3n functionals are computed at 
each iteration level by the present method. Also, to compute a Jacobian 
(nxn) , n^ elements must be stored and a total number of N arithmetic 
operations (addition, subtraction, multiplication) must be done where 
N = n![2 + 1/2! + 1/3! + ... + l/(n-l)!]. In comparison with this the 
present method computes and stores n-diagonal elements of the Jacobian. 

There are certain interpolatory methods available to solve stiff 0DE*s 
[9]. These algorithms, developed by Certaine and Jain reduce differential 
equations into integral equations which were integrated by approximating 
the integrands by interpolation polynomials over the range of integration. 
Although both of these methods have high accuracy and are A-s table, for 
large systems they are not practical. 
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The present method, having a simple algorithm, is applicable to most 
nonlinear systems in general. It can be applied to nonlinear integral equa- 
tions or integro-dif ferential equations. Some applications with regard to 
this are given in [12] . 

5. Applications to Nonlinear Stiff ODE^s 

Ex; 1 X = 50/x - 50x , x(0) = 

The analytical solution is x(t) = / 1 + exp (-100 O. To solve this 

stiff ODE we approximate x by a backward difference formula 

X. = X, , + At(50/x. - 50x.) 

J J-1 33 


If we replace x. 


by 


X and X. 

J-1 


by xO we get 


X = At(50/x - 50x) + xO 

Then, g(x) = At(50/x “ 50x) + xO 

Also, 

g(x^) - g(x*) = -50 At/—— + 

\xV 

where x = the value of the root and x^ = value of x at some kth 
iteration. 

Writing g(x^) - g(x*) = ~ x*) , we note that a^ is a D-element 

— 3 

if At < 10 . However, if we define a new 

g(x) = X + a( (At(50/x - 50 x) + xO) - x) 



Then, 


a^= (1 - a) - 50 


a At f — + l) 


and a, is a D-element if a 
k 


is sufficiently small and positive. 


-1 


Computationally, it has been found that for a = 1.0 and At = 10 the 

— 6 3 

method failed, whereas for a = 10 and At = 10^ it did not. 
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Ex; 2 X = -10004 x + 10000 y** 

y = -y + X - 
x(0) = y(0) = 1 

This system of stiff ODE's is given in [3] . It has an approximate 
analytical solution given by 


y = 


10004 e~^*^ 
10008 - 4 e"^^ 


nV3 


10000 4 

^ =Ioo^y 


Using backward Euler's difference scheme, this nonlinear system may be 
expressed as: 


X. = At(-10004 X. + 10000 y'^) + x. , 
J J J-1 


yj - 






X 

O 


= y. 


= 1 . 


Replacing x^ by x, y^ by y, x^_^ by xO and x^_^ by yO, 
we get: 

X = At(-10004x + 10000 y**) + xO 


y = At(-y + X - y^*) + yO 

Let us express the system as x = F(x,y) and y = G(x,y); where 

F(x,y) = X + {At(-10004 x + 10000 y**) + xO - x} 

G(x,y) = y + 02 {At(-y + x - y^*) + yO - y} 


At some (k+l) iteration level: 

k+l * ,/k *\.,,k 

X - X = bii(x - X ) + bi2(y - y ) 

k+l * , k+l *\ , , y k 

y - y = a2i(x - X*) + b22(y - y ) 


where 


bii = 1 - Oj - 10004 At Oj 
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bi 2 = lOOOOAtOjCy^ + y*^)(y + y*) 


321 = At • 02 


b 22 = 1 - 02 - At 02 (y^ + y*^) (y + y*) 


From (19) , 

0 o“ 


bn bi2 






321 0 


CM 

CM 

o ^ 

Here 

I 

^11 

bl2 


\ - V‘ \ - 


-321 bn 


-321 bi2 + b22 


Case: 1 Transient solutions 

Choosing = a 2 = 0.1 and At = 0.005 results obtained by using per- 
turbed functionals have been presented in table 5. 

Case : 2 Steady-State solution 

If the objective is to compute solutions at the steady state, variable 
time steps may be used and At could be very large. Using a sequence of 
At = 1, 10, 10^, 10^... 10® it has been found that after a total of 55 
iterations steady state solutions were found using Radio Shack TRS-80 color 
microcomputer (16K). Computational time was about one minute with x*s 
converging to zero faster than y’s. Newton’s method did not show this mode 
of convergence. 

Chemical Kinetic Models . 

Magee and Chat ter jee [4] developed models on chemical kinetics and 
sought solutions of these problems which should be time-accurate. These 
models are nonlinear and consist of a sequence of stiff ODE’s. Some of 
these chemical species grow from a concentration of 10 moles/liter at 
t = 10 ^ sec to 10 ^ moles/liter at t = 1.0 sec; whereas by the same time 
certain other species which have values of the order 10”^® moles/liter at 
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t = 10 sec grow up to 10 moles/liter at t = 1.0 and stay almost 
unchanged. Since the equations are stiff, in the transient process numerical 
effects of some components of the solution which decay much faster in com- 
parison with other components are quite difficult to capture in a computa- 
tional process. This leads to some inaccuracies in the solution which 
eventually causes instabilities. For example, if the concentration of a 
species is found to be -10 mole/liter but whose true concentration is 
10 moles/liter, this causes severe instabilities. In such computations 
which are very ’^sensitive” with regard to small errors, principles of per- 
turbations could be applied in order to stabilize the algorithm. Indeed 
this was the finding when two distinct chemical kinetics problems were solved 
by the method. Let us consider them, and discuss some of the very interest- 
ing computational findings verified by actual experiments. 

Models representing irradiation of water with y-rays have been consid- 
ered. Two conditions for water were taken (i) Acid Water and (ii) Pure Water. 
Ex 3 : Irradiation of Acid Water 

The y-radiation energy is absorbed with the creation of molecular 
(H 2 , H 2 O 2 ) and radical (H, OH) products which may be treated as if formed 
homogeneously in the system. The mechanisms related to formations of these 
species are described in [5,6]. Seven species created by the radiation 
participate in thermal reactions, summarized in Table 1. If the irradiation 
is continuous at the rate of 100 I’s electron volts per liter per sec, the 
differential equations which describe the concentration changes are given 
in Table 2. These equations are stiff . Here creation terms give concen- 
trations of created species in terms of moles per liter per second. 

The system described in Table 2 may be expressed as: 
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where x = [(H) (OH) (H 2 O) (H 2 ) (H 2 O 2 ) (HO 2 ) ( 02 )]^ and f = [fj £2 ... 

where fi corresponds to the right side of the first equation, f 2 corre- 
sponds to the right side of the second equation, etc. Applying a backward 
Euler-dif f erence scheme we get 

x*^^^ - x^ - hf(x^^^) =0, h = At 

This nonlinear system for may now be expressed as: 

F(X) = 0 

where X = x'^'*'^ which could be put in the form as 

X = G (X) 
o 

where = X + a • F(X) and a = diag (a^, 02 , ... ay). As discussed 

in the Ex: 1, 2, 0 < < 1. 

A numerical solution is found for continuous irradiation at constant 
rate which starts with the molecular concentrations of H 2 , H 2 O 2 and O 2 equal 
to zero. For this case the radical concentrations approach "stationary 
values" very quickly and as known from the chemical nature of the system 
that H^ and are destroyed in a chain reaction, these concentrations 

which build up linearly at first, approach stationary values on a longer 
time scale. Computational results showed these behaviors of various concen- 
trations of this complicated chemical process. Figures 1, 2, 3 show these 
results for various values of I. (More results on this project will be 
published from Lawrence Berkeley Lab.) From "zero" values, H and OH quickly 

“ft „Q 

approached 2 x 10 and 2.5 x 10 , respectively, in 0.2 sec and then slowly 

decreased attaining "stationary values" at time equal to 1 sec. However, 

H 2 , H 2 O 2 , HO 2 and O 2 attain their stationary values after almost 15 sec. 


14 



Chatter jee and Magee [4] applied an implicit second-order Runge-Kutta 
method to solve this problem and recorded its failure. 

Ex; 4 Irradiation of Pure Water 

The treatment of pure water is somewhat more complicated simply because 
it involves more species. Here is so low in concentration that the 
hydrated electrons are not converted into H atoms before the track reactions 
occur and charged species must be treated explicitly. There are eleven 
equations giving the rates of change of concentrations of eleven species. 

In Table 3, the thermal reactions and in Table 4 the differential equations 
of the system are given. Continuous irradiation starting with zero decomposi- 
tions for all concentrations of species excepting H 2 O = 55 moles/liter, 

H 30 ^ = l.OE - 07 moles/liter and OH = l.OE - 07 moles/liter was done. At 
this stage, the ODE's were represented as integral equations which were 
approximated by trapezoidal rule. The method virtually failed when the 
derivatives were approximated by backward Euler's difference formula. In 
the code a very stringent convergence criterion, namely, 

WMAX = max |a)^| < 10 

3 

3 

was used. Since this should give a very high order of accuracy, it was felt 

that such a criterion is necessary in order to get an almost perfect "mass 

balance" and "charge balance" both of which should be (theoretically) equal 

to "zero" at all times. In submicroseconds, the radical species like H, OH, 

®Aq rapidly than the molecular species O 2 , O 2 , H 2 , H 2 O 2 . 

These have strong impact upon the mathematical model for the subsequent 

—8 

chemical yields. Thus At = 10 was chosen initially. It was noticed that 
while all other species were growing, O 2 and O 2 stayed zero up to 5 x 10~® sec. 
Then O 2 grew faster and O 2 > O 2 at all time levels. These computational 
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properties of these species which possibly have some strong impact upon the 
solution as time increases, cannot be detected if a large time step is used. 

In Fig. 4 we see the transient stages of the growths of some of the 
species up to 10 ^ sec. Up to this stage 1130"^ and OH ^ did not show any 
changes from their initial concentrations. However in Fig. 5 some notice- 
able changes of their values were found around t = 1.0 sec. The values of 
these concentrations seem to be in agreement with those known both theoreti- 
cally and experimentally. 

In Fig. 6 we see all the species reaching steady-state around 
t = 30 secs. 

For a further checking of the validity of these computational results 
an additional computer-run was taken with an initial value of H2O2 = 10 
These should cause a faster growth of O2, one of the most sensitive species, 
and its steady-state-concentration should be identically the same as that of 
H30'^. This vital aspect of the chemical reaction process is clearly depicted 
in Fig. 7. Also it was noticed that whereas O2 should increase H2O2 must 
decrease and both should merge at a steady-state configuration. This result, 
as evidenced by experiment, was also found computationally. 

6. Discussions 

k+1 k 

From equation (12) it is clear that if 9^0^ ’ = 1, the method fails. 
It has been proved and demonstrated in [14] , that the perturbation parameters 
Stabilize the algorithm and speed up the rate of convergence provided the 
Jacobian of the matrix representation of the system has nonzero diagonal 
elements. When diagonal terms are null, artifically they are brought in and 
are damped out as convergence is approached. This may be called nonlinear 
scaling. However, if diagonal elements are zeros and nonlinear scaling is 
not performed the rate of convergence is significantly decreased and the 
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algorithm could even be destabilized. This was found with regard to solu- 
tion of several problems. The essential strength of the method lies in 
having functionals with diagonal nonlinearity. Fortunately, the set of 
equations given in Table 2 and Table 4 fulfill these requirements. 

By computer experimentation it has been found that for nonlinear sys- 
tems having multiple roots, the present scheme is not quite effective •* 

Although on a trial basis this method can be applied to solve a non- 
linear system, it may be useful to apply the properties of D-Mappings 
to the functionals in order that the method may be applied successfully 

[8,15] . Such an analysis is generally quite complicated. 

7. Concluding Remarks 

This is a preliminary report. The objective with regard to the solu- 
tion of the chemical kinetic problem was primarily to study the treatment 
growths of the species, not the final equilibrium values which were some- 
what known to the chemists [4] through experimental data. With this regard, 
the present method served the purpose of the researchers very well. 
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Table 1: Reactions in Irradiated Water (Acid) 



Reactions 

Reaction rate 
constant, 
£/(mol s) 


A. Recombination of Primary Radicals 

1. 

H + H ^ H2 

1 X 10^0 

2. 

H + OH -> H 2 O 

2.4 X iQlO 

3. 

^ OH + OH -»• H 2 O 2 

4 X 10^ 


B. Reactions of Radicals 

and Product Molecules 

4. 

H + H202^H20 + OH 

1 X 10^ 

5. 

OH + H 2 O 2 -> H 2 O + HO 2 

5 X 10^ 

6. 

OH + H 2 H 2 O + H 

6 X 10^ 

7. 

HO 2 + H H 2 O 2 

1 X 10^° 

8. 

HO 2 + OH H 2 O + O 2 

0 

I — 1 

0 

tH 

X 

9. 

HO 2 + HO 2 ^ H 2 O 2 + O 2 

2 X 10^ 

10. 

H + O 2 ^ HO 2 

1 X lO^O 
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Table 2: Differential Equations for Transient Species 

and Radiation Products in Irradiated Water (Acid) 


^ (H) = 3.711 - 2ki(H)2 - k2(H)(0H) - k4(H)(H202) + kg (OH) (H 2 ) - k7(H02)(H) - kio(H)(02) 

^ (OH) = 2.951 - k2(H)(0H) - 2k3(0H)2 + k4(H)(H202) - k5(0H)(H202) - k6(0H)(H2) - kg(H02)(0H) 

^ (H 2 O) = -4.511 + k2(H)(0H) + k4(H)(H202) + ks(0H)(H202) + ke(0H)(H2) + k0(HO2)(OH) 

^ (H 2 ) = 0.401 + ki(H)2 - ke(0H)(H2) 

^ (H 2 O 2 ) = 0.781 + ka(0H)2 - k4(H)(H202) - k5(0H)(H202) + k7(H02)(H) + k9(H02)^ 

^ (HO 2 ) = ks(0H)(H202) - k7(H02)(H) - ks(H02)(0H) - 2kg(H02)2 + kio(H)(02) 

^ (O 2 ) = k8(H02)(0H) + k9(H02)2 - kio(H)(02) 


The range of integration is from 0 sec to 20 sec. 



Table 3: Reactions in Irradiated Neutral Water 


Reaction rate 
constant 

Reactions il/(inol s) 


A. Recombination of Primary Radicals 


1. 

H + H H2 

2. 

e". + H Ho + OH" 

Aq 

3 . 

e". + e“, -s- Ho + 20 H" 

Aq Aq 

4 . 

e". + OH ^ OH" 

Aq 

5. 

H + OH -> H2O 

6. 

OH + OH -V H2O2 

7. 

HaO'*' + H 

8. 

H3S + OH- H2O 


B. Reactions of Radicals 

9. 

H + H2O2 -»■ H2O + OH 

10. 

+ H2O2 -> OH + OH" 

11. 

OH + H2O2 H2O + HO2 

12. 

OH + H2 H2O + H 

13 . 

HO2 + H H2O2 

14 . 

e", +02"^ O2" 

Aq 


1 X IQlO 
2.5 X IQlO 
6 X 1Q9 

3 X lO^O 
2.4 X IqIO 

4 X 10^ 

2.3 X lO^O 
3 X 10^0 

with Product Molecules 

1 X 1Q8 

1.2 X IqIO 

5 X 10^ 

6 X 10^ 

1 X IQlO 
1.9 X IQlO 
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Table 3: Reactions in Irradiated Neutral Water 

(cont.) 


Reaction rate 
constant 

Reactions £/ (mol s) 


B. Reactions of Radicals with Product Molecules 

(cont.) 


15. 

HO2 + OH H2O + O2 

1 

X 

10^0 

16. 

HO2 + HO2 H2O2 + O2 

2 

X 

106 

17. 

H + O2 -> HO2 

1 

X 

0 

0 

00 
I — 1 

02 ~ + H^O HO2 + H2O 

3 

X 

0 

> — 1 

0 

T— 1 


C. Dissociation Reactions 

+ 

19. H 2 O H 3 O + OH 

+ 

20. HO 2 -> H 3 O + 02 “ 

*Rate constant, sec~^ 


5.5 X 10-6* 
1 X 106* 


23 




Table 4: Differential Equations for Transient Species and 

Radiation Products in Irradiated Neutral Water 


d(H) 

dt 


= - 2ki(H)2 -k2(e"^q)(H) - k5(H)(0H) + k7 (nto) (e"^^) - kg(H)(H202) + ki2(0H)(H2> 

- ki3(H02)(H) - ki7(H)(02) 


^ ■■ - - H(e-^,)(OH) - k, (HsS) - k, „ (Oj) 

^ (OH) = - k4(e-^^)(0H) - k5(H)(0H) - 2k6(0H)2 + kg(H)(H202) + kj q ( e“^^) (H 2 O 2 ) - kj 1 (OH) (H 2 O 2 ) 

- ki2(0H)(H2) - kis(H02)(0H) 

^ (H^O) = - k7(H30)(e"^^) - k8(Hto)(OH") - kj g (Hg^) (O 2 ") + kig(H20) + k2o(H02) 

(H 2 O) = ks(H)(OH) + kg(H3S)(OH“) + kg(H)(H202) + ki 1 (OH) (H 2 O 2 ) + ki2(0H)(H2) + kis(H02)(OH) 

+ kig(H3$)(02 ) - kig(H20) 

(H 2 ) = ki(H)2 + k2(e"^^)(H) + k3(e"^^)2 - ki2(0H)(H2) 

^ (H 2 O 2 ) = kg (OH) 2 - kg(H)(H202) - ki g (e"^^) (H 2 O 2 ) - ki 1 (OH) (H 2 O 2 ) + ki3(H02)(H) + kig(H02)2 


+ 0,55 I 
+ 2.65 I 

+ 2.70 I 
+ 2.65 I 

- 4.10 I 
+ 0.45 I 


+ 0.70 I 


Table 4: Differential Equations for Transient Species and 

Radiation Products in Irradiated Neutral Water 

(cont .) 


^ (oh") = k2(e~^^)(H) + 2kg (e“^^) 2 + k4(e“^^)(0H) - keCHgOHOH-) + ki o (e“^^) (H2O2) + kig(H20) 

(HO2) = kii(0H)(H202) - ki3(H02)(H) - ki5(H02)(0H) - 2ki6(H02)2 + ki7(H)(02) + kj 3 (H3O) (O2") 

^ (O2) = - ki4(e-^^)(02) + ki5(H02)(0H) + ki6(H02)2 - ki7(H)(02) 

^ (02") = ki4(e"^^)(02) - ki8(H30)(02") + k2o(H02> 


- k2o(H02> 



Table 5 


At = 

0.005, k = No. 

of iterations 

for convergence 

at a given time 

step. 

t 

Exact Solution 

Dey 

k 

X 

y 

X 

y 

0 

1.0 

1.0 

1.0 

1.0 

0 

0.25 

0.368021198 

0.778953673 

0.368799110 

0.779286476 

5 

0.50 

0.135369435 

0.606629568 

0.136012224 

0.607287244 

4 

0.75 

0.049796505 

0.472436166 

0.050160955 

0.473250544 

3 

1.00 

0.0183185671 

0.367930928 

0.018499225 

0.368797576 

3 

1.25 

0.73174825E-03 

0.286467519 

6.82246374E-03 

0.287398821 

3 

1.50 

2.47909732E-03 

0.223160237 

2.516106E-03 

0.223965900 

3 

1.75 

9.12006087E-04 

0.173797232 

9.27933E-04 

0.174533507 

2 

2.0 

3.35507793E-04 

0.13535337 

3.42219E-04 

0.136011530 

2 

2.25 

1.23426334E-04 

0.105413292 

1.26209E-04 

0.105991891 

2 

2.5 

4.54059955E-05 

0.0820959478 

4.6546E-05 

0.082598004 

2 


26 


Appendix: A 

A flow chart to solve: X = F(X,Y), Y = G(X,Y) using perturbed functionals. 
Notations : = Partial Derivative of F with respect to X 

Gy = Partial Derivative of G with respect to Y 
KMAX = Maximun No. of iterations 
EP = Epsilon 

WX = Perturbation applied to X 
WY = Perturbation applied to Y 
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Figure Captions 


1. Figure 1. Concentrations of Species (Ex: 3) vs. time for 
I = 6.667 X 10 ^ in the logarithmic scale. 

2. Figure 2. Concentrations of Species (Ex: 3) vs. time for 

_7 

I = 6.667 X 10 in the logarithmic scale. 

3. Figure 3. Concentrations of Species (Ex: 3) vs. time for 

_8 

I = 6.667 X 10 in the logarithmic scale. 

4. Figure 4. Concentrations of Species (Ex: 4) vs. time for 

_7 

I = 6.667 X 10 in the logarithmic scale, up to t = 3 x 10 

5. Figure 5. Concentrations of Species (Ex: 4) vs. time for 

—7 

I = 6.667 X 10 in the logarithmic scale up to t = 3 secs. 

6. Figure 6. Concentrations of Species (Ex: 4) vs. time for 

-7 

I = 6.667 X 10 in the logarithmic scale up to t = 

(Here steady state is reached for all the species.) 

7. Figure 7. Concentrations of Species (Ex: 4) vs. time 

—7 

I = 6.667 X 10 in the logarithmic scale up to t = 

(Initial value of H 2 O 2 = 10 ^ moles/liter.) 


30 secs, 
for 

30 secs. 


sec. 
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CONCENTRATION (moles/l) 



TIME (sec) 


XBL 823-8495 









CONCENTRATION (moles/I) 



Fig. 3 
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MOLES/LITER 



TIME(sec) 


XBL 824-9063 


Fig. 7 
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